#' Estimate the percentage of green duckweed within a red circle 
#' 
#' @aliases leaf_areas leaf_area
#' 
#' @usage 
#' leaf_area("the image file name including the path to file")
#' leaf_areas("the fold containing all the images only")
#' 
#' @param image the file name of the image
#' @param fold the fold that contains jpg image only
#' @param plotfinal whether to save the final identified figure to a file with the same name
#' 
#' @details 
#' This function is designed to calculate the percentage of green area within the red circle.
#' Performance of the function critically depends on the following things:
#' 
#' 1. take the photo as clear as possible
#' 
#' 2. make the color of the circle as red and isotropic as possible
#' 
#' 3. remove other alga, especially green alga, as far as your can
#' 
#' 4. reflection of red and green color by the glass beaker will introduce bias too
#' 
#' @note 
#' The jpeg package is required to read the jpg/jpeg images. The spatstat pacakge is required to plot the final figures.
#' 
#' @return 
#' The percentage of green area in the red circle
#' 


leaf_areas=function(fold,plotfinal=TRUE){
  images=dir(fold)
  p=unlist(lapply(images,leaf_area,plotfinal=plotfinal))
  return(p)
}

leaf_area=function(image,plotfinal=TRUE){
  
  require("jpeg")
  
  pc=readJPEG(image)
  
  #remove black
  black=pc[,,1]<0.4 & pc[,,2]<0.4 & pc[,,3]<0.4
  #remove white
  white=pc[,,1]>0.7 & pc[,,2]>0.7 & pc[,,3]>0.7
  
  #remove small RGB diff areas
  other=abs(pc[,,1]-pc[,,2])<0.1 & abs(pc[,,3]-pc[,,2])<0.1
  
  
  noti= (black | white | other)
  pc[,,1][noti]=NA
  pc[,,2][noti]=NA
  pc[,,3][noti]=NA
  
  #circle= (pc[,,2] <0.3)  
  circle=abs(pc[,,3]-pc[,,2])<0.1
  circle[is.na(circle)]=0
  circle=circle+0
  
  yq=quantile(which(apply(circle,1,sum)>10))
  xq=quantile(which(apply(circle,2,sum)>10))
  
  xsr=xq[c(2,4)]
  ysr=yq[c(2,4)]
  
  nnc=which(circle==1)
  xc=nnc%%dim(circle)[1]
  xc[xc==0]=dim(circle)[1]
  yc=(nnc-xc)/dim(circle)[1]+1
  
  loc=data.frame(x=yc,y=xc)
  test=data.frame()
  j=1
  testcircle=function(n)
  {
    i=n
    while(j<=i)
    {
      selrow=sample(row(loc),3)
      x1=loc[selrow[1],]$x
      y1=loc[selrow[1],]$y
      x2=loc[selrow[2],]$x
      y2=loc[selrow[2],]$y
      x3=loc[selrow[3],]$x
      y3=loc[selrow[3],]$y
      a=2*(x2-x1)
      b=2*(y2-y1)
      c=x2^2+y2^2-x1^2-y1^2
      d=2*(x3-x2)
      e=2*(y3-y2)
      f=x3^2+y3^2-x2^2-y2^2
      l1=sqrt((x1-x2)^2+(y1-y2)^2)
      l2=sqrt((x2-x3)^2+(y2-y3)^2)
      l3=sqrt((x1-x3)^2+(y1-y3)^2)
      if(l1>50&l2>50&l3>50)
      {
        suml=l1+l2+l3
        xpot=(c*e-b*f)/(a*e-b*d)
        ypot=(c*d-a*f)/(b*d-a*e)
        r=sqrt((xpot-x1)^2+(ypot-y1)^2)
        test1=data.frame(xpot=xpot,ypot=ypot,ra=r,suml=suml)
        test=rbind(test,test1)
        j=j+1
      }
    }
    return(test)
  }
  
  #plot(im(circle))
  testre=testcircle(1000)
  #points(x=testre$xpot,y=testre$ypot,pch=".")
  xdensity=density(testre$xpot,from=xsr[1],to=xsr[2])
  ydensity=density(testre$ypot,from=ysr[1],to=ysr[2])
  xp=xdensity$x[which(xdensity$y==max(xdensity$y))]
  yp=ydensity$x[which(ydensity$y==max(ydensity$y))]
  #points(xp,yp,col=2,pch=19)
  
  dd=sqrt((testre$xpot-xp)^2+(testre$ypot-yp)^2)
  rr=testre[which.min(dd),]$ra
  #abline(h=yp-rr)
  #abline(h=yp+rr)
  #abline(v=xp+rr)
  #abline(v=xp-rr)
  #TODO: adjust the rr to 
  
  
  
  #remove the most irrelavent areas
  subpc=pc[round(yp-rr):round(yp+rr),round(xp-rr):round(xp+rr),]
  subcircle=circle[round(yp-rr):round(yp+rr),round(xp-rr):round(xp+rr)]
  dimyx=dim(subcircle)
  focusyx=dimyx/2
  newx=rep(1:dimyx[2],each=dimyx[1])
  newy=rep(1:dimyx[1],time=dimyx[2])
  
  dd2=sqrt((newx-focusyx[2])^2+(newy-focusyx[1])^2)
  
  #try to estimate a more accurate radius by small reducing dd2 in each step
  
  
  
  outcircle=dd2>rr
  dim(outcircle)=c(dimyx[2],dimyx[1])
  outcircle=t(outcircle)
  subpc[,,1][outcircle | subcircle]=NA
  subpc[,,2][outcircle | subcircle]=NA
  subpc[,,3][outcircle | subcircle]=NA
  
  #by color
  #duck=(subpc[,,1]<0.6 & subpc[,,3]<0.3)+0
  duck=abs(subpc[,,2]-subpc[,,3])>0.2+0
  duck[duck!=1]=NA
  #duck=(subpc[,,1]>0.2 & subpc[,,1]<0.6 & subpc[,,3]>0 & subpc[,,3]<0.3)+0
  
  #finally gives the percentage of duck coverage in the container.
  re=sum(duck,na.rm=T)/sum(!(outcircle | subcircle))
  
  if(plotfinal){
  require(spatstat)
  jpeg(filename = paste(image,".jpg",sep=""))
  plot(im(subcircle))
  plot(im(duck),add=T)
  dev.off()
  
  }
  return(re)
}
